
# delim ;  
clear all;
set matsize 11000 ;
eststo clear ;
capture program drop add_lab_FE;
capture log close ; 
log using "../logs/optimality_conditions.log", text replace;

*settings some settings that are common for the set of plots on this dofile; 

local pdf_plot_settings = 	`"
				xmtick(1992(1)2008)
				xlabel(1992(4)2008)
				ylabel(#3,)
				ymticks(##5)
				scale(*1.5) 
				xtitle("")
				legend(off)"';
*scale(*1.3): This rescales the font size of the whole figure;	
********************************************************************;



	

*parameters ;
local iri_max=180;
global delta =0.1;
global delta10 =0.1;

global kappa=0.0000016;


*------------------------------------------------------------------------------;
*1. Get Parameters for the parametric approaximation;
*------------------------------------------------------------------------------;
use "../data/hpms_hwy_stats_88_08.dta", clear  ;


*---------------------------------------------------;
* - Generate variables and general cleanings ;
*---------------------------------------------------;
do "_programs_sample";
do "_setup_expenditure.do";
do "_cleaning_and_new_variables_sample.do";

	gen I2_t_exp_I2Xtime=I2_t_exp_I2*time;
	gen I2_t_app_I2Xtime=I2_t_app_I*time;
	label var I2_t_exp_I2Xtime "$ \mathds{1}_{ist}(q)y^Q_{st}\times t$";

** generate a variable of average number of years between maintenance in a segment;

bysort id (year): egen N_maint=total(I2_t);
bysort id (year): egen xx=count(year);
gen av_year_maint=xx/(N_maint+1);


label var av_year_maint "mean t for I in seg. i";
tab av_year_maint;
drop  xx;

eststo clear;
	eststo, prefix(table_) : reghdfe  Diri  c.I2_t_exp_I2 time  I2_t_exp_I2Xtime   , abs(id ) cluster(`cluster_var') ;
	sum_stats;
	
local A1=_b[I2_t_exp_I2];
local A2=_b[I2_t_exp_I2Xtime];

tab year time;
*------------------------------------------------;
*II. CONSTRUCTION;
*------------------------------------------------;
*Define price index ;
clear;
use "../data/state_year_all_80_15.dta", clear  ;

* - Generate variables and general cleanings ;
*---------------------------------------------------;

do "_setup_expenditure.do";
do "_cleaning_and_new_variables.do";
do "_programs_universe";

tab year periods ;

drop if exp_L_IH_SF12a_r==.;
drop if exp_L_IH_SF12a_r==0;
drop if exp_L_IH_SF12a_r<0;
drop time2;
gen time2=time^2;
label var time2 "\( t^2\)";
gen exp_L_IH_SF12a_rXtime=exp_L_IH_SF12a_r*time;
gen exp_L_IH_SF12a_rXtime2=exp_L_IH_SF12a_r*time2;


gen time_92=time-8 if year>1991;
tab time_92 year;

gen exp_L_IH_SF12a_rXtime_92=exp_L_IH_SF12a_r*time_92;
local IV L4FHWA_apport_LW_r;
gen `IV'Xtime =`IV'*time;
gen `IV'Xtime_92 =`IV'*time_92;
gen `IV'_sample=1 if `IV'!=.;
sum app_IH_r exp_L_IH_SF12a_r;

*eststo ,prefix(filter_): reg Du_lane_miles_IH 	 time   exp_L_IH_SF12a_r exp_L_IH_SF12a_rXtime  io36.state  if year>1991,   cluster(state);
eststo ,prefix(filter_): 	eststo ,prefix(filter_): reg Du_lane_miles_IH 	 time_92   exp_L_IH_SF12a_r exp_L_IH_SF12a_rXtime_92 if year>1991 ,   cluster(state);
local Acons1=_b[exp_L_IH_SF12a_r];
local Acons2=_b[exp_L_IH_SF12a_rXtime];


eststo ,prefix(filter_): 	eststo ,prefix(filter_): reg Du_lane_miles_IH 	 time time2  exp_L_IH_SF12a_r exp_L_IH_SF12a_rXtime  exp_L_IH_SF12a_rXtime2  ,   cluster(state);
local Acons1_all_sq=_b[exp_L_IH_SF12a_r];
local Acons2_all_sq=_b[exp_L_IH_SF12a_rXtime];
local Acons3_all_sq=_b[exp_L_IH_SF12a_rXtime2];

eststo ,prefix(filter_): ivreg2  Du_lane_miles_IH 	 time   (exp_L_IH_SF12a_r exp_L_IH_SF12a_rXtime=`IV' `IV'Xtime),   robust first  ffirst  saverf saverfprefix(rf_3) savefirst  savefprefix(fs_3) cluster(state);

local Acons1_IV_all=_b[exp_L_IH_SF12a_r];
local Acons2_IV_all=_b[exp_L_IH_SF12a_rXtime];

eststo ,prefix(filter_): ivreg2  Du_lane_miles_IH 	 time_92   (exp_L_IH_SF12a_r exp_L_IH_SF12a_rXtime_92=`IV' `IV'Xtime_92) if year>1991,   robust first  ffirst  saverf saverfprefix(rf_3) savefirst  savefprefix(fs_3) cluster(state);
local Acons1_IV_92_08=_b[exp_L_IH_SF12a_r];
local Acons2_IV_92_08=_b[exp_L_IH_SF12a_rXtime];



eststo ,prefix(filter_): ivreg2  Du_lane_miles_IH 	 time_92   (exp_L_IH_SF12a_r exp_L_IH_SF12a_rXtime_92=`IV' `IV'Xtime_92) if year>1991,   robust first  ffirst  saverf saverfprefix(rf_3) savefirst  savefprefix(fs_3) cluster(state);
local Acons1_IV_92_08=_b[exp_L_IH_SF12a_r];
local Acons2_IV_92_08=_b[exp_L_IH_SF12a_rXtime];




*------------------------------------------------------------------------------;
*2. Calculate optimality conditions
*------------------------------------------------------------------------------;



use "../intermediate_data/national_vars" , clear ;
tset year ;

foreach var in  rev_grand_total_mt_r rev_total_h_r rev_total_mt_r	rev_total_m_fuel_h_r	{;	   
gen `var'_IH=	`var'*(s_vmt_IH 	/s_vmt_NHS)  ;
gen `var'_IH_per_VMT=(`var'_IH*10^6)/s_vmt_IH   ;
*NOTE WE MULTIPLE BY 10^6 to get the revenue in dollars and not thousands USD ;
};

gen tcons=[_n]-1 ;
gen tcons2=tcons^2;
gen t=[_n]-12 if year>1991;
tab t year;
gen P_const_92_08=`Acons1'+`Acons2'*t;
gen P_const_all_sq=`Acons1_all_sq'+`Acons2_all_sq'*tcons +`Acons3_all_sq'*tcons2;
gen P_const_IV_all=`Acons1_IV_all'+`Acons2_IV_all'*tcons;
gen P_const_IV_92_08=`Acons1_IV_92_08'+`Acons2_IV_92_08'*t;
gen P_resurf=`A1'+`A2'*t;




gen q_t =s_iri_IHnat if year>1991 ;
label var q_t "Interstate quality (inches)";
*ren u_vmt_IH u_vmt_IH;
ren u_lane_miles_IH L_t;
*gen  r=cond(DGS10_r>0,DGS10_r/100,.);
gen r_raw=DGS10_r/100;
reg r_raw year , ;
predict r, xb;

# d ;
tw (scatter  P_const_beta year if year>1983,msize(*0.5) ) (lowess  P_const_beta year if year>1983 ,msize(*0.5) lw(dot) lw(*1.3) color(gs0))
(line P_const_92_08 P_const_IV_all P_const_IV_92_08 year  if year>1983,msize(*0.5) lp(dash solid dot) lw(*1.8 *1.1 *3 ) )
,ytitle("") xtitle("") 
`pdf_plot_settings' 
leg(on) leg(on) 
xlabel(1984(4)2008)
xmtick(1984(1)2008)
leg(region(color(none)) col(1) ring(0) pos(2)
		size(*0.85)
		
label(1 "Non parametric")
label(2 "Non parametric (Smooth)")
label(3 "Baseline")
label( 4 "IV All ")
label( 5 "IV 92-08 ")		
		
		)
 ;
graph export "${output}/figures//FigureB5_P_const.pdf", replace ;


lowess  P_const_beta year ,msize(*0.5) generate(P_const_lowess) ;


gen r_h=rev_total_m_fuel_h_r_IH_per_VMT; 
gen p_L_92_08=1/P_const_92_08*10^6;
gen p_L_IV_all=1/P_const_IV_all*10^6;
gen p_L_IV_92_08=1/P_const_IV_92_08*10^6;
gen p_L_all_sq=1/P_const_all_sq*10^6;
gen p_L_beta=1/P_const_beta*10^6;
gen p_L_lowess=1/P_const_lowess*10^6;

gen p_q_beta=-1/P_resurf_beta*10^6;          
gen p_q=-1/P_resurf*10^6;
gen p_m=1/P_maint_beta;
gen m_h=exp_IH_mtn_SF12a_r*10^6/u_vmt_IH;


# d ;
tw (scatter  p_L_beta year if year>1983 ,msize(*0.5) ) (lowess p_L_beta year if year>1983,msize(*0.5) lw(dot) lw(*1.3) color(gs0) )
(line p_L_92_08 p_L_IV_all p_L_IV_92_08 year if year>1983 ,msize(*0.5) msize(*0.5) lp(dash solid dot) lw(*1.8 *1.1 *3 ) )
,ytitle("") xtitle("") 
`pdf_plot_settings' 
xlabel(1984(4)2008)
xmtick(1984(1)2008)
leg(region(color(none)) col(1) ring(0) pos(10)
		size(*0.85)		
label(1 "Non parametric")
label(2 "Non parametric (Smooth)")
label(3 "Baseline")
label( 4 "IV All ")
label( 5 "IV 92-08 ")	)
ylabel(0(10000000)50000000)
 ;
graph export "${output}/figures//FigureB5_p_L.pdf", replace ;


egen m_h_mean=mean(m_h);



gen LHS_mtn=r_h;
gen LHS_basic=r_h;

sum p_q if year==1993;
gen p_q93  = r(mean) if year>1991 ;

tw (scatter  r_ra year if year>1991,msize(*0.5) )
(lowess  r_ra year if year>1991,msize(*0.5) )
(line r  year if year>1991,msize(*0.5) )
,ytitle("") xtitle("") 
`pdf_plot_settings' 
 ;
graph export "${output}/figures//FigureB4_i_rate.pdf", replace ;



tw (con LHS_mtn year if year>1991,msize(*0.5) ),ytitle("") xtitle("") 
`pdf_plot_settings' 
ylabel(0(0.01)0.02) ;
graph export "${output}/figures//Figure6_opt_mtn_LHS.pdf", replace ;


*Generate variables ; 
*-------------------------------; 
foreach x in lowess all_sq beta 92_08 IV_92_08 IV_all {;
    local p_L p_L_`x';

gen RHS_basic_`x'=(r*`p_L'*L_t + r*p_q*q_t*L_t+(${kappa}*p_q*u_vmt_IH))/u_vmt_IH;

gen RHS_mtn_`x'=  m_h + RHS_basic_`x';

sum RHS_mtn_`x'  , d;
gen ratio_basic_`x'=LHS_basic/RHS_basic_`x';
gen gap_basic_`x'=LHS_basic-RHS_basic_`x';


gen ratio_mtn_`x'=LHS_mtn/RHS_mtn_`x';
gen gap_mtn_`x'=LHS_mtn-RHS_mtn_`x';



	 

	*------------------------------------------;	
	* Evaluating the Euler Equations; 
	*------------------------------------------;
	sort year ;

	*- Create prices to create counterfactual;
	*----------------------------------------;



	sum `p_L' if year==1993;
	gen `p_L'93  = r(mean) if year>1991 ;





	foreach a in 1 4{;

		global alpha=0.`a';

		* - Euler condition on L:
		*---------------------------;
		* individual parts; 
		gen r_h_L0_alpha`a'_`x'= L_t/((1-${alpha})*u_vmt_IH);

		gen r_h_L1_alpha`a'_`x'   =((1+r[_n-1])*`p_L'[_n-1]-`p_L'       + p_q*(q_t + ${kappa}*(1-${alpha})*(u_vmt_IH/L_t) -q_t[_n+1]));
		*i) pq93;
		gen r_h_L1_pq93_alpha`a'_`x'   =((1+r[_n-1])*`p_L'[_n-1]-`p_L'       + p_q93*(q_t + ${kappa}*(1-${alpha})*(u_vmt_IH/L_t) -q_t[_n+1]));

		*ii) pl93;
		gen r_h_L1_pl93_alpha`a'_`x'   =((1+r[_n-1])*`p_L'93[_n-1]-`p_L'93       + p_q*(q_t + ${kappa}*(1-${alpha})*(u_vmt_IH/L_t) -q_t[_n+1]));

		*-Euler equation;
		gen r_h_L_alpha`a'_`x'= m_h +r_h_L0_alpha`a'_`x'*r_h_L1_alpha`a'_`x';
		
	

		* - Euler condition on q:
		*---------------------------;
		gen r_h_q0_alpha`a'_`x'= q_t/(${alpha}*u_vmt_IH);

		gen r_h_q1_alpha`a'_`x'   =(p_q[_n-1]*L_t[_n-1]*(1+r[_n-1])      -p_q*L_t*(1-${kappa}*${alpha}*(u_vmt_IH/(q_t*L_t))));
		
		gen r_h_q_alpha`a'_`x'=m_h +r_h_q0_alpha`a'_`x' *r_h_q1_alpha`a'_`x';	
		
		


	};
	

	

  ;


};

*Figures;

		tw (con  RHS_mtn_92_08 year if year>1991,msize(*0.5) ),  
		ytitle("") xtitle("") `pdf_plot_settings' ylabel(,format(%3.2fc)) ;
		graph export "${output}/figures//Figure5_opt_mtn_RHS_92_08.pdf", replace ;
	
			tw (con r_h_q_alpha1_92_08 r_h_q_alpha4_92_08 year if year>1991, msize(*0.5  *0.5) lp( dot solid ) color(gs0%40 gs0%90 ) ),
		legend(label(1 "{&alpha}= 0.1") label(2 "{&alpha}=0.4") 

		region(color(none)) col(2) ring(0) pos(7)
		size(*0.8))	ytitle("") xtitle("") `pdf_plot_settings'
		xtitle()
		ylabel(,format(%3.2fc))
		 leg(on) ;
		graph export "${output}/figures//FigureB6_r_h_q_alphas_92_08.pdf", replace ;


		tw (con r_h_L_alpha1_92_08 r_h_L_alpha4_92_08 year if year>1991, msize(*0.5  *0.5) lp(dot solid)   color(gs0%40 gs0%90 ) ),
		legend(label(1 "{&alpha}= 0.1") label(2 "{&alpha}=0.4") 	
		region(color(none)) col(2) ring(0) pos(7)
		size(*0.8))	ytitle("") xtitle("") `pdf_plot_settings'
		ylabel(,format(%3.2fc))
		 leg(on) 
		;
		graph export "${output}/figures//FigureB6_r_h_L_alphas_all_92_08.pdf", replace ;



# d ;

foreach y in  p_q p_L_92_08 u_vmt_IH r{;
	sum `y' if year==1992;
	gen `y'_92  = r(mean) if year>1991;
};

gen r_high=r+0.04;

gen RHS_basic_92_08_check =(r*p_L_92_08*L_t + r*p_q*q_t*L_t+(${kappa}*p_q*u_vmt_IH))/u_vmt_IH;
gen RHS_basic_92_08_r_high =(r_high*p_L_92_08*L_t + r_high*p_q*q_t*L_t+(${kappa}*p_q*u_vmt_IH))/u_vmt_IH;
gen RHS_basic_92_08_r92 =(r_92*p_L_92_08*L_t + r_92*p_q*q_t*L_t+(${kappa}*p_q*u_vmt_IH))/u_vmt_IH;
gen RHS_basic_92_08_p_q92 =(r*p_L_92_08*L_t + r*p_q_92*q_t*L_t+(${kappa}*p_q_92*u_vmt_IH))/u_vmt_IH;
gen RHS_basic_92_08_p_L92 =(r*p_L_92_08_92*L_t + r*p_q*q_t*L_t+(${kappa}*p_q*u_vmt_IH))/u_vmt_IH;
gen RHS_basic_92_08_VMT92=(r*p_L_92_08*L_t + r*p_q*q_t*L_t+(${kappa}*p_q*u_vmt_IH_92))/u_vmt_IH_92;

sum RHS_basic_92_08_check RHS_basic_92_08;
pause;

foreach c in r_high r92 p_q92 p_L92 VMT92{;

	gen RHS_mtn_92_08_`c'=  m_h + RHS_basic_92_08_`c';



	 #  d ;
	 pause ;
};


*------------------------------------------;
*Counterfactuals
*------------------------------------------;
# d ;
preserve;

	keep if year==1992 | year==2007;
	
	global var_euler_cond   r_h_q_alpha1_lowess r_h_q_alpha1_beta r_h_q_alpha1_92_08 r_h_q_alpha1_IV_92_08 r_h_q_alpha1_IV_all 
							r_h_q_alpha4_lowess r_h_q_alpha4_beta r_h_q_alpha4_92_08 r_h_q_alpha4_IV_92_08 r_h_q_alpha4_IV_all 	
							r_h_L_alpha1_lowess r_h_L_alpha1_beta r_h_L_alpha1_92_08 r_h_L_alpha1_IV_92_08 r_h_L_alpha1_IV_all 
							r_h_L_alpha4_lowess  r_h_L_alpha4_beta r_h_L_alpha4_92_08 r_h_L_alpha4_IV_92_08 r_h_L_alpha4_IV_all ;
					
	global var_euler RHS_mtn_lowess RHS_mtn_beta RHS_mtn_92_08 RHS_mtn_IV_92_08 RHS_mtn_IV_all;
	
	global var_counterfactual RHS_mtn_92_08_r_high RHS_mtn_92_08_r92  RHS_mtn_92_08_p_q92  RHS_mtn_92_08_p_L92  RHS_mtn_92_08_VMT92;
	
	keep year ${var_euler} ${var_counterfactual}
	;
	
	gen id=1;
	reshape wide ${var_euler} ${var_counterfactual} , j(year) i(id);
	
	ren RHS_mtn_92_08_* Counterfactual_*;
	ren *1992 y1992* ;
	ren *2007 y2007* ;
	

	reshape long y2007 y1992,j(var) i(id) string;
	drop id;


	
	gen ratio=y2007/y1992;

	gen eq="A_Counterfactual" if 
								var=="Counterfactual_r92" | 
								var=="Counterfactual_p_q92" | 
								var=="Counterfactual_p_L92" |
								var=="Counterfactual_VMT92";
								
	replace eq="B_Sensitivity" if var=="RHS_mtn_IV_all" |
								   var=="RHS_mtn_IV_92_08" |
								   var=="RHS_mtn_beta" |
								   var=="RHS_mtn_lowess"
								   
	;
	
	replace eq="__"  if var=="RHS_mtn_92_08" | var=="Counterfactual_r_high" ;
	
	replace var="AA_RHS_mtn_92_08" if var=="RHS_mtn_92_08";
	replace var="AA_RHS_mtn_92_08_r_high" if var=="Counterfactual_r_high";
	pause;
	
	sort var;

	 mkmat y2007 y1992 ratio ,matrix(table_CF) rownames(var) roweq(eq); 
pause;
 # d ;
 
	*change tex ->  csv if you want to see on excel;
	esttab matrix(table_CF, fmt(%9.2fc %9.2fc  %9.2fc )) 
	using "${output}/tables/Table5.tex", replace 
	
		substitute(
		"&    table_CF&            &            \\" "" 
		ratio "2007/1992"
		y1992 "1992"
		y2007 "2007"
		RHS_mtn_lowess   "Non parametric (Smooth)"
		RHS_mtn_beta     "Non parametric"
		RHS_mtn_IV_92_08 "\hline IV 92-08"
		RHS_mtn_IV_all   "IV All"
		Counterfactual_r92 "\(r_{92}\)"
		Counterfactual_p_q92 "\(p^q_{92}\)"
		Counterfactual_p_L92 "\(p^L_{92}\)"
		Counterfactual_VMT92 "\hline \(\sc{VMT}_{92}\)"
		AA_RHS_mtn_92_08       "Baseline"
		AA_RHS_mtn_92_08_r_high       "Baseline (r high)"
		"A_" "A. "
		"B_" "B. "
		"__          &            &            &            \\" ""
			
			)

				prehead(\begin{table}[!htb]\centering
\caption{Sensitivity and Counterfactuals \label{tab:Counterfactuals}}
\begin{tabular}{l*{4}{c}}
\hline\hline
)

postfoot(\hline\hline
\end{tabular}
\flushleft
\noteTableCounterfactuals
\end{table});
				
				
				
	
	
	
	
restore;



*------------------------------------------;	
* Construction figure adding the line used in the 
Calibaration Exercise. 
--> This figure is in the paper; 
*------------------------------------------;
merge 1:1 year using  "../intermediate_data/construction_betas_1.dta"; 
*This table is created in construction_regressions.; 

drop if year<1984;

twoway 
 (rcap up5 down5 year, lstyle(ci))
 (scatter beta year, m(D))

 (line P_const_92_08 year, lp(solid) lw(*1.5) color(gs0%70))

 ,
 
 `pdf_plot_settings' xtitle("") 
  xlabel(1984(4)2008)
  xmtick(1984(1)2008)
  ; 
 graph export "${output}/figures/Figure3_yearXexpL_0_linear.pdf", replace ;
 
 ren up up_const ;
 ren down down_const;
 ren beta beta_const;
 drop _merge;
 *------------------------------------------;	
* Construction figure adding the line used in the 
Calibaration Exercise. 
--> This figure is in the paper; 
*------------------------------------------;
# d;
merge 1:1 year using  "../intermediate_data/resurf_betas_1.dta"; 
*This table is created in construction_regressions.; 

drop if year<1984;

twoway 
 (rcap up5 down5 year if year>1991, lstyle(ci))
 (scatter beta year if year>1991, m(D))

 (line P_resurf year if year>1991, lp(solid) lw(*1.5) color(gs0%70))

 ,
 
 `pdf_plot_settings' xtitle("") 
  ; 
 graph export "${output}/figures/Figure2_year_X_I2exp_1_linear.pdf", replace ;

*------------------------------------------;	
* Table With the time series; 
*------------------------------------------;

preserve ;
	keep year u_vmt_IH L_t q_t r_h r p_L_92_08 p_q m_h ;

	format u_vmt_IH L_t q_t  p_L   %9.2fc  ;
	format L_t %11.1fc;
	format m_h r_h r %9.3fc;
	format p_q %9.1fc;

	replace u_vmt_IH=u_vmt_IH/10^9;
	replace p_L=p_L/10^6;

	 mkmat  u_vmt_IH L_t q_t r_h r p_L p_q m_h ,matrix(table) rownames(year); 

	 
		*change tex ->  csv if you want to see on excel;
		esttab matrix(table, fmt(%9.2fc %9.1fc %9.2fc %9.3fc  %9.3fc  %9.2fc %9.1fc %9.3fc)) 
		using "${output}/tables/TableB5_national_variables.tex", replace 
			substitute(\_ _ 
						"u_vmt_IH" "VMT"
						"&       table&            &            &            &            &            &            &            \\" ""
									"&    VMT&         L_t&         q_t&         r_h&           r&   p_L_92_08&         p_q&         m_h\\"
				"&    \(\textsc{vmt} \times 10^9\)&         \(L_t\)&        \(q_t\)&         \(r_h\)&           \(r\)&        \(p^L \times 10^6\)&         \(p^q\)&         \(m_h\)\\")

					prehead(\begin{table}[!htb]\centering
	\caption{National variables for the calibration \label{tab:national_vars}}
	\begin{tabular}{l*{8}{c}}
	\hline\hline
	)

	postfoot(\hline\hline
	\end{tabular}
	\flushleft
	\noteTableVars
	\end{table})
					
					
					
					
					
				;
 
restore;
	
